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Abstract 

In the paper the flow in a thin tubular structure is considered. The velocity of the flow stands for 
a coefficient in the diffusion-convection equation set in the thin structure. An asymptotic expansion 
of solution is constructed. This expansion is used further for justification of an asymptotic domain 
decomposition strategy essentially reducing the memory and the time of the code. A numerical 
solution obtained by this strategy is compared to the numerical solution obtained by a direct FEM 
computation. 
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1 Introduction 

The paper is devoted to the strategy of numerical implementation of the asymptotic partial decomposition 
of the domain for the tubular structures of a complicated geometry. We will consider the Stokes flow 
in this structure and the convection-diffusion and sorption process for some diluted substance. First we 
consider an asymptotic expansion of the solution. We emphasize the importance of the boundary layers 
in the neighborhood of some special structural elements of the tubular domain, such as the bifurcations 
of canals and " stenosis areas" . That is why some multiscale strategy should be applied to the analysis 
of the convection-diffusion process: the ID limit description in the canals will be coupled with some 2D 
zooms in these special structural elements. 

In section 2 we define a tubular structure as a union of thin rectangles connected by some domains 
of small diameters. The Stokes equation and the diffusion-convection equation are set in this domain. 
For the Stokes equation the Dirichlet conditions are respected at the lateral boundary with some given 
inflow and outflow. For the diffusion-convection equation we pose the Robin type condition at the 



lateral boundary with some given inflow and outflow concentrations. The viscosity and the diffusion are 
constant out of some "stenosis area" where they may have variations. The varying viscosity can be used 
for the modeling of a clot in the blood circulation process. Indeed, if at some part of the domain the 
viscosity is great, then, applying the idea of the fictitious domain method, we can exclude this part of 
the domain from the flow i'll'Oi'l (see Remark [lj. 

In section 3 we consider the Stokes equation in tubular structure. The asymptotic expansion of the 
solution for constant viscosity has been obtained in [T] . In this section we construct the boundary layer 
correctors for the varying viscosity in the stenosis areas. 

In section 4 the convection-diffusion equation is considered. First we construct the asymptotic ex- 
pansion in an infinite tube (subsection 4.1). Then using this expansion as a regular ansatz we add the 
boundary layer correctors in the stenosis zones (subsection 4.2) in a bifurcation area and in the en- 
trance/exit elements (subsections 4.3 and 4.4). In the subsection 4.5 the leading term of the asymptotic 
expansion is presented. The justification of the asymptotic expansion follows the scheme: estimate for 
the residuals and application of the a priori estimates for the initial problem. 

Section 5 describes one version of the partial asymptotic domain decomposition strategy for the mass 
transport problem in a tubular structure. 

Finally, section 6 develops the numerical experiment comparing the direct numerical solution of 
the 2D problem and the asymptotic solution of the partially decomposed problem. The results of this 
experiment confirm good coincidence of the exact solution and the approximate solution obtained by the 
method of asymptotic partial decomposition of domain. 



2 Geometry of tubular structure and setting of the problem 

We will introduce the tubular domain which consists of three types of structural elements: canals, 
bifurcations and stenosis areas. This tubular structure is similar to the rod structures introduced in 
[2] and the tube structures or pipe structures introduced in [3]; we consider a new element that is, the 
stenosis area, simulated by varying coefficients of the equation (viscosity and diffusion coefficients) and 
not by geometric singularity. 

Let us remind the definition of a tube structure. 

Let ei, e„ be n closed segments in R 2 which have a single common point (i.e. the origin of the 
coordinate system) and let it be the common end point of all these segments. Let 9i,...,6 n € (0,1) 
be n positive numbers. Making a change of variables (rotation) such that the new axis x\ denoted x e { 
contains the segment and the second new axis x% % is orthogonal to e*, we define 

B\ = f (an.xa) e M 2 : xf € (0, |e*|) , x% G 





Let 70 be a bounded domain containing O, ji be a bounded domain containing Oi, the end point of a 
(different from O). We assume for simplicity that diam(7^) , diam(7o) < 2. Consider the homothctic 
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contraction of 70 in — times with the center of the homothety in O and denote Jq the image of 70 by 
this homothety, i.e. 



75 = {(x 1 , 2 ; 2 )eM 2 : 6 7o} 



e e 

In the same way we consider 



{x 1 ,x 2 ) - Oi 



7? = l^r,^) G K : + Oi e 7< 

the homothetic contraction of 7$ in - times with the center of the homothety in Oi. 
Define the one bundle tube structure 

Be = (ysf) u ( u H 

such that it is a connected domain with the C 2 — smooth boundary. 

In a more general case, we will consider several one-bundle structures B e i, B £ 2, B £m such that 
every of these structures is associated to some segments: 

en, ei ni for B £l , 
esi, e2« 2 for £? £2 , 

^mlj ■••) &mn n for B £nl . 

Assume that if two of these segments have common point then it is an end point for each of these 

m 

segments; assume that the union (J i? e? is a connected domain with C 2 — smooth boundary. In this case 

9=1 

m 

we will call this union B e = {J B eq a multi-bundle tube structure. If an end point of some segment 

9 =1 

is not an end point for all other segments then such end point will be called solitary. 

We will consider the Stokes equation and the convection-diffusion equation in such tube structure. 
The boundary condition is the vanishing velocity for the Stokes equation everywhere except of some 
special parts of the boundary (entrance and exit). These parts are some connected parts Ti of the 
boundary of smoothing domains 7j. We assume that the end points of segments corresponding to this 7^ 
are solitary end points. Let Ti, ...,r r be these parts of the boundary. 

We consider the Stokes equation in such a tube structure with the varying viscosity coefficient /1 : 

div (x) (Vu £ + (Vu 6 ) T ) - p s /) = U 0) , 



div \ 



(1) 



where the divergence is taken with respect to the elements of each line of the matrix \i e (x) (vu £ + (Vit e 
p £ I, and the convection-diffusion equation 

- div (K £ (x) Vc e ) + u £ (x) ■ Vc e = g (xf ) . (2) 

Assume that g = in some neighborhood of the end points of the segments, g £ C k+2 (e;) for all 
segments ej. 

We will assume that fi £ and K £ arc positive constants fi and h respectively, everywhere except some 
"stenosis areas" where they have a form: 

H £ (x) = n + M ( X - Xs ) and K £ (x) = x + K [ % - Xs ) (3) 
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where M and K arc measurable bounded function having a finite support inside the ball B 2 J with 

the center — and the radius 2, such that. 

e 

3xi > : n £ (x) , K e (x) > x\. 

Here x B are some points belonging to the segments of the graph of the structure, they are different 
from the end points and are independent of s. 

The sorption will be modeled by the boundary condition for the diffusion-convection equation, i.e. 
let us consider the boundary conditions: 

u e = on the lateral boundary dB £ \ ( |J T t J , (4) 

dc e ( r \ 

K s (x) = e(3c E on the lateral boundary dB e \ ( U r t J , (5) 

u e = G(^^j onr t (6) 
c e = qt = const on T*, t = 1, r (7) 
where Xb t is an end point, inside 7t, of a corresponding segment, n is an outer normal, G € Cq (7 t ) and 



Remark 1 The varying viscosity and diffusion coefficients can be used for the modelling of a clot in the 
blood circulation process. Let us remind the fictitious domain method. Consider an example: the Poisson 
equation Au = f posed in a bounded domain G with the boundary condition u\qg = 0; / € L 2 (G) . The 
fictitious domain method reduces this problem to the problem set in a larger rectangular Rd G : 

div (K u (x) Vu w ) = F (x) , x S R, 

Uu\dR = 



where 



1, x e G, . v _ f f(x) , x G G, 

" ^'-^ u; ^ x G R\G, [ { ~ \ 0, x G i?\G. 




Fig. 2 



For smooth dG, one can prove that, as lu — > +oo 

0, a; G R\G, mH {R) " 

5o £/ie yieW ti^ vanishes in the fictitious part of the domain. The same effect holds for the Stokes equation, 
where the varying viscosity can be used to modelling the absence of flow in the fictitious part of domain 
occupied by the clot. 



4 



3 The Stokes equation 



We will consider separately the problem for the Stokes equation and the convection-diffusion one. First 
we apply the results of [3] and get the asymptotic expansion of the solution. At the second stage we 
assume that the velocity u e is known and consider the convection-diffusion equation with the velocity 
coefficient corresponding to the first term of the asymptotic approximation. 

The asymptotic solution of the Stokes problem was considered in pQ. The only difference is related 
to the "stenosis areas" where the boundary layers are constructed as follows. 

A stenosis area can be simulated by a varying viscosity in some close neighborhood of the origin 
of the coordinate system. Then we can consider one channel parallel to the Ox\ axis, i.e. G £ = 

(— f , f ) x ^— ^, . We will not take care of the ends of this channel because these ends are supposed to 
be junction points with some other channels and the construction of the bifurcation boundary layers is 
described in [T] . So we will try to construct a solution of problem (JTJ) stabilizing to a Poiscuillc solution 

as > +oo. 

£ ^ 

Let equations {T]) be considered in this channel G e and let the viscosity coefficient \x e have a structure 

where /i > is constant and M has a support inside the ball B (0, 2) = {£ G R 2 : £ 2 + £ 2 < 4} such that 
M (£) > 0. 

Then if the right-hand side f e is equal to zero, the asymptotic solution out of the boundary layer 
zone (at some finite distance from 0) is a Poiscuillc flow: 

«p = c 1 ^( 2;2 o _T )^, p(x)=c lXl +c 2 , (8) 

where c±, C2 are some constants. 

Then the boundary layer corrector has a form (e 2 U (£) , eP (£)) and (U,P) is a solution of the 
following problem 

div 4 (( M + M (f )) (v c tf + (V C C/) T ) - P (0 /) = -ci div c (m (0 ( i I ) 6 
div 5 {7 = 0, £ G (—oo, +oo) x 

where the right hand side 

- Cl div 5 (m(£)( ! g 

has a support inside B (0, 2) . It is well known ([H [5]) that this problem has a unique solution ({/, P) 
stabilizing to (0, const) at the infinity. Then (u p + e 2 U (-) ,p(x) + eP (-)) satisfies equation (JTJ) and 
flow coincides with the Poiseuille flow at a finite distance from zero with an exponentially small error 
O (e~^) , a > 0. If P (£) — ► for £i — > — oo and P (£) — > c + for £i — * +oo, then we have to proceed the 
following special gluing of pressure at zero: 
I) we redefine a new (discontinuous) 

c\x\ + C2, for xi < 0, 
c\X\ + c 2 + ec + , for xi > 0; 



I I 

2' 2 



p(x) -- 

2) we redefine a new (discontinuous) (denoted by the same letter): 



P(^.-l P (0, for6<0, 
1 j ' I P(0-<M-. for £i >0. 
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Mention that the sum p (x) + eP (f ) is still smooth and this new P (£) — > for £ x — > ±00, i.e. it has a 
standard boundary layer shape. The justification of this expansion follows laterally [1] 

4 The convection-diffusion equation 

Consider the diffusion-convection problem |(3J), and in a tube structure B e , where the coefficient 
K £ is given by formula Q and u £ in ([2]) is replaced by the given vector-valued function V £ having the 
following structure: 



Vs(x) 



Vi(2=fL) fovx€B s : \s=ft\ <2, i = 0,..,n, 
V. (s=Sm.) for x e B e : MN < 2, 



( Vp \ ' J) for xeB e :x e 2 > £ e h \^=f i \ > 2, |^^| > 2, for all i and s, 



where = — ((£2*) — 4) an d (0 an d V* (£) are some given smooth vector- valued functions 



with finite support in the ball B (0,2) , and x s are the "stenosis nodes" and Oi are the ends of the 
segments e,. 

Really, the structure of the velocity field u e is more complicated: out of the balls of radius 2e 
surrounding the nodes Oi and the stenosis points x s the velocity differs from the Poiseuille flow by some 
exponentially decaying boundary layer functions. Here we simplify the structure of the velocity field 
replacing u e by V e . So we consider here the problem: 

- div (K e (1) Vc e ) + V c (x) ■ Vc £ = g {x e {) (9) 

with the boundary conditions 

K e (x) — e(3c E on the lateral boundary dB £ \ ( U ^tj 1 (10) 
c e — It = const on T t , t = 1, r. (11) 



dc 

We assume that p is a constant and /3 < 0. It means that the sorption takes place: the outflow — jt 



c V 



dn 

is positive. It is well known that if at least one of given concentrations c s is equal to zero, then the 
Poincarc-Fricdrichs inequality for c s holds with a constant independent of e. Then it could be shown that 
there exists a constant uq, independent of e, such that, if ||u £ |j ioo ^ B ^ < uq then there exists a unique so- 
lution c e of problem (fT2|) — (|T3 |) . It satisfies an a priori estimate \\c e \\^(b s ) — c {\\9\\l 2 (b e ) + ll^lli? 1 (B s 
where c is independent of e (see [6]). 

4.1 The asymptotic expansion in a channel 

Consider the equation 

- div (xVc Er ) + V (—J ■ Vc er = 3 (zi) , (12) 

in the infinite channel G e = (-00, +00) x (-§, §) ,where V (&) = f ^ ^ ^ , y p (£2) = - (£§ - |) ; 
we consider the boundary conditions 

x—^=e(3c £r if 2:2 =±-. (13) 
on 2 

We will construct a function 

k 

cto(xi,b) = ^e?c j {x u &) (14) 



G 



such that it satisfies ([12]). (fT3|) up to the terms of order e k if g £ C fc+2 (M). Substituting ([14]) into (fl~2l) . 
(fH| we get 

_ 2 / 9 2 c,- (xi,£ 2 ) . 9 2 c 7 _ 2 T , , t . 9c,_ 2 



where i2g is a discrepancy 



j=k-l,k 

and if £2 = ±5, 



2 • 

j=0 ^ 2 " j=0 

Equating the terms of the same power of e, we get 

9 2 c 7 - 9 2 C 7 _ 2 TA ,„ . 9c 7 _2 / v r 

*b% = ~*^r + p (6) + 5 632 

and for £ 2 = ±^ 

9c 

±X^| 6=± i =0Cj-a (iCi.fij). 

The necessary and sufficient condition of existence of c,- : the condition 

5 9 2 c 7 ^ 9c,- / 1\ 9c, / 1 
-^72~«S2 = x— \xi, - - x^— I an, -- 
1 9^ 96 V 2 / 9 6 V 2 

implies 

■ , /„, u .dCj- 2 \ , •, p ./ / 1 



9. 2 ' \ V * (6) -feT/ + 5 ( ' Tl) Sjl = * ^' 2 J + Cj - 2 ^ ' - 

1 

where (■) = -d£ 2 - 

Each function c,- is sought as a sum Cj (an, £2) = c,- (xi) + c,- (xi,^) , (c,-) = and (c,-) = c,-. So for 
Cj_2 we get the equation: 

So, we get an algorithm of the successive determination of c,- , Cj : 

9% 



3 = 



" ^2 -°> & G (-5, 1) 
x^| 6=±4 = 0| (co} = 



and so, cq = 0; 

-^^ + (^)£+5(^i) = 2/3co; (15) 
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j > 1; find c,- : 



'96 

Remark ci = 0, and find Cj 



x-^ = -^- + V p ( 6 )-^- + fl (x 1 )5 J1> 6c(-U) (i6) 

±X ^/I«2=±* = ^ C J'-2 fcl'SO , (Cj) = 0. 



5s 2 Ox\ \ p d Xl / ' J ' ' V V 2/ V 2 

Lemma 2 i?ac/i Cj (xi,^) *s a polynomial function in £2 of degree 2j. 

Proof. By induction we prove that Cj is a polynomial function of the degree 2j. Indeed, V p is a 
quadratic function and so the right hand side (fT6|) is a polynomial of order 2 (j — 2) + 2. And so, after 
two integrations of (fT6|) we check that is a polynomial of order 2 (j — 2) + 2 + 2 = 2j. ■ 



So we have constructed an asymptotic approximation Q14p which satisfies equation (I12p up to the 
remainder and conditions ((T3|) up to the remainder e fc /? (^1,^2) |^ 2 =±i- 

Let us eliminate this remainder in the boundary conditions. To this end we will add a corrector 
e k+1 Ck+i (^1,^2) such that 



dx 2 

i.e. 



k+1 c k+ i (xi, y) = e k+2 /3c k +i (xi, y) + £k P c k (xi, — ) for x 2 = ±|, 



±^--e fc B fe+1 (.T 1 ,6) =e fc+2 /3c fe+1 (x 1 ,6) + £ fc /3c fc ( a;i ,6) for & = (17) 



For example, 



* 9 ,- \ ,- \ 9 / \ ,- 

i\ f . 1 V „ / i\ 1 \ i\ „ / 1 



a+i = j (6 - 2 J (6 + 2 J - ) - ( 

satisfies (fTT)l . Then cir — £ fc+1 Q..+i denoted by ci£ satisfies the boundary conditions (fT3j) exactly and 
the equation (|T2"]) up to the remainder + R^ , where 



R [k) = l-div (xVc fc+1 ) - V p (6) ^} 



4.2 Structural element "stenosis area". 

Consider the convection-diffusion equation in the channel G £ with the modified coefficients x and V : 
they are replaced by the functions 

K e (x) = x + K (~) and V e (x) = + V (~) 

such that J\ e (a;) > >fi > for all x, divV (-) = and the support of the functions K (£) and V(£) 
belongs to the ball B (O, 2) = e I 2 : |£| < 2} . 
Consider the convection-diffusion equation 

-dw(X e (x)Vc B )+V e (^).Vc e =ff e (xi) > inG e (18) 
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with the boundary conditions 

K E (x)-^=ef3c E7 x 2 = ±-. (19) 

Let us construct a function c|g satisfying equation (|18[) and conditions Q19p up to remainder of order 
e k . To this end we will consider the asymptotic approximation ci^ of the previous section completed by 

(k) 

the special boundary layer corrector c esbl , having the form 



such that 



k 

„(*) . 
"esbl 

j=0 



|C/,-(OI<cie- c ^l. 

Denote 

L £ = div (K £ V) - V e ■ V. 

Define 

f 1 for |t|>2; 
PW \ for |t| < 1; 

such that p and p' are bounded. Consider the following asymptotic approximation 

k 



Then 



LJS - 9 = -L. (eg) (l - P (y))) - .9 + (l - P (y)) 



where — p (y — ^j \ Lec£r — an d |-Rfe| < e k const. 



Expand cir in powers of e£i : 



Let us denote K (£) the function x + A' (£) and V (£) the function Co) +V (£) . We have 

^ ( x ^ p (t))) = E £J " 2 ^ (* (0 v € (E f ^gr 1 (0.6) >< (! - p(&» 



E ^ v (0 • ^ (E I ^ (o. 6) >< (i - p 



J=0 



So, for Uj we get equation 

div 6 (K (0 VsUj (0) - v (0 ■ V e ^_i (0 + T, (0 =0, 6 6 R, 6 € 



1 1 

2' 2 
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where 



Ti (0 = div 6 (k (£) V 4 ( £ |^ (0, 6) x (1 - p &)) j 



For the boundary conditions, in the same way 
fe ^ rr k 



± ^'-^ (o H ± f; e*-^ (o |g o, ^ - 1) = efi ( E ^ + E ^ o» (&) - 1) ) . & = 4 



3=0 

and so 
where 

Sj (0 



0$ 



J'=0 



±£*(0 



^i 9 ^ (0,6) (i - p(&)) - Ifr (0,6 



Necessary and sufficient condition of existence of a bounded solution C/j : 

{Tj (0 ~ V (0 • V s ^_ x } = {^} + + {S j }_ 

where 



)(i- P (6))-/3C/ J - 



{■} 



and 



<-w. 



(_oo,0)x(-i |) 



(-oo,0) 



(0,+oo)x(-i i) 

■I?2 = ±i^l- 

(O. + oo) 



This condition gives one interface condition for Q J 1 (0 + ) and Q J 1 (0~) 



dx i 



^ (0^2 



9xi 



(0-)- / A" (Odf: 



9xi 



(0+)- 



fe ( o-6)^-/_>(e)^i(o +J 6)^ + 



+ all other terms of Tj except Z = 1 for div. 
So, if we pose Qj = f_\ # (£) (0~, 6) «a - /ji # (0 &i (0+, 6) dfc, we get condition: 



9xi 



^(0,6)^6 



9xi 



<Jj 



(20) 



where 



fJj 



{Tj (0 - v (f) ■ v,r. .! - Qj + {Sj}, + {S j }_ , 



where in Tj (£) there are not the term corresponding to I = 1 for dw. 
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Now we solve the problem on Uj exponentially stabilizing to some constants. We choose this constant 
equal to at — oo. At +00 we have Uj — > q~j. To make Uj — > as £1 — > +00, we subtract this constant 
for all £1 > 0; we set 

f Uj for £i < 
3 \ Uj- Qj for ^ > 0. 

This function Uj is exponentially decaying at 00. Then we put the compensating condition for Cj at 
xi = : [cj] = <&•. 

Another interface condition for Cj is (|20jl . 

4.3 The boundary layer in a "bifurcation area". 

Consider the convection-diffusion equation in a one bundle structure with the common point O (i.e. the 
origin of the coordinate system). We will construct an asymptotic expansion in the neighborhood of the 
point O (that will be good in all the channels Bf). First we construct a regular expansion as in section 
14.11 for every branch Bf making the local change of variables by rotation in such a way that new x\ axis 
(denoted x^) contains the segment ej. Let us denote the approximation c£p for ej as 

w-±* (*$). 

Here {x\\ x 1 ^') are new local variables related to the segment e,;, and is the thickness of the channel 
Bf. Then, as in section we construct the boundary layer corrector 

3=0 

such that 

\Uj(0\ <cie-«l«l, ci,c 2 > 0. 
Consider an asymptotic approximation 

i— 1 ^ ^ \ / j— 

where Xi is a characteristic function of Bf (xi ( x ) = 1 if x G and if not); a is a number such that 
each point of the circle of a radius ae with center O belongs to at most one rectangle Bf (i = 1, ...,n) 
and it does not intersect the domain "/q . It means that the variation of the function p from 1 to occurs 
out of this circle. 




Fig. 3 
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Assume that the right hand side is equal to zero in some neighborhood of O. Then for every i = 1, n 
for x G Bf 



L £ C eb 



= -l. ^ (i - , + (i - , " ') + 



k 

' X 



where i?fe = p \ ^ 1 J ^Leci? 12 ' — 3 J • Inside the circle B (0;ae) we g 

Here the terms ^1 — p ^ aXl ^ ^L 6 ci? e * — 5^ and p ^ aXl ^ ^L e ci? e ' — <7^ are bounded by e k const. 
So we will define C/j in such a way that out of the circle -B (0; ae) we have 



u =o ° / »=i 



to be equal to zero up to the terms of order e k and inside the circle, 

k 

X 

e 



ME^(J)I=°- 



Expand each cl? 6 ' in powers of 



9 (a^ 



fc r f l d l c ei = T e 

= E £r ElirrS( <^)+^, €?- 



e 



Consider domain fi = 70 U \J fij, where 0, is the half-strip |(£i>£2) € R 2 : ^ > 0, |£f*| < y j- . So for 
?7j we get equation: 



xAt^ (£) - v (o ■ v e f/ 3 -_i (0 + Tj (0 = 0, e 6 n 0l 

where 

T?* (£) if £ 6 Hi out of the circle B (0, a) , 
inside the circle B (0, a) , 

- (£*)' 3V* 



Tj (0 



Tf ft) = xAt ( E ^ ^ (0,ff) >< (1 -P«'D 



,1=0 9(x e {) 

-v(e)-v,(E^^^(o,^) 
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For the boundary conditions on 8Qq: 



: E" ^ ■ E ' - !) 



dc* 



3=0 



3=0 



and 



k k 

= eP^e^U-j + ^e j cf (p (a^) - 1) , out of the circle B (0, a) 

3=0 3=0 
k QJJ . k 

± e^~ x >c—^- = e/3 inside the circle B (0, a) 

3=o n t 3=0 



up to the terms of order e*. This gives 



dU 



where 



^ (0 = 



f 3 r^V d l+1 c £z 



(O^j'Kl^fO) Jf/ r 2, S out of the circle B (0, a) on ffl„ 



+/3Uj-2, inside the circle B (0, a) on dilo 
Necessary and sufficient condition of existence of a bounded solution Uj : 

{2i(0-V(0-V € U J -_i} no ={5 i } ano) 

where 

{•}a = / < and {-}aa = / - ds - 

dc ei . dc ei 

This condition gives one interface condition for (0 ei )i the limit value of the derivative on ei 

in the origin of e,; : 

n ^e, n ggH 

{Tj} n = - V x6»j — ^ (0 e! ) - XI TT^" (° £! > ) + a11 other tcrms of t j except 1 = 1 for A « ■ 

i=i ^i* »=i dx i* 



So, we get Kirchoff type condition: 



where gj depends on c , ...,c,_2. 

Now we solve the problem on Uj in Qq; it stabilizes exponentially to some constants qji at every 
branch f2^. This solution Uj corresponds to some values of c^' (0 ei ) that enter in the expression for 
Tj' (when I — 0). If we change these values adding qji then the solution Uj will be transformed into 
Uj — p (atii*) qji which tend to zero as |£| — > +oo. 

It means that for every j wc have to solve first the problem for Cj, that is, an ordinary differential 

equation on every e.j, Kirchoff condition for 1 in and continuity condition for c 7 in 0, solve the 

ox i 

problem for Uj and then modify the values Cj (0 £i ) in such a way that Uj — ► as |£| — » +oo. 
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4.4 The entrance/exit element 

In this case the boundary layer is constructed in the following way. It satisfies equation (|12[) . boundary 
conditions (| 13[) everywhere except the part T s , where Uj = q s 8oj- There is no more necessary and suffi- 
cient condition for the existence of a bounded solution. This bounded solution exponentially stabilizes 
to some constant. Then we redefine the value Cj (xe, s ) in such a way that U — ► as — > +00. In 
particular, c (x bs ) = q s and U Q = (1 - p) g s . 



4.5 Algorithm of assembling of an asymptotic solution 

Let us describe more precisely the algorithm of the assembling of the different structural elements in the 
construction of the asymptotic expansion. 

For co we get from the construction of the subsection 14.11 the ID convection diffusion equation 



(y P )e,^ + 9 i9(xi)=2f^, for every e,, (22) 



At every stenosis point x s we get the interface conditions 

= and [cq 1 ] = 0. (23) 



rfr ei 



dx{' 

At every bifurcation point Xb that is an end point of segments e\, e„ we get 

X E^W = ( 24 ) 

z— 1 1 

and Cq ! (x;,) = Cq 3 (xb) Vi, j 6 {1, n} (continuity). 

At every entrance/exit point xt we set c^* (xt) = qt, t = 1, .., r, where art is an end-point corresponding 

to r t . 

Then for every j we get analogous problem for Cj with some right hand sides depending on previous 
approximations co, ...,cy_i. We solve the problem for Cj and then for Uj. We define the interface values 
Cj* (xb) for all segments e, having as an end point in such a way that Uj — ► for |£| — > +00. 

Remark 3 Every function Cq 1 related to e, depends on the local variable x & { such that x e { = in one 
of the end points of e^. Therefore we get two local variables "starting" from each of two ends of the 
segment a; therefore, if we denote x\ i+ and x e {~ these two different local variables, we get the relation 
for these variables x\ i+ = |e, | — x\ l ~ . Let us give more details. Every segment having two end points 
Ojj and Oi 2 can be associated to two possible local coordinate systems: one of them ( denoted by Ox ei i ) 
has its origin in Oi 1 and another in Oi 2 (denoted by Ox £i z ). Therefore for every differential equation 
fSjjJ) or for every junction condition this system should be chosen and fixed. The evident change variable 
relation for any function f defined on the segment e% is: 

f(x^)=f{\e i \-x e 1 i *), 

where |e,| is the length of the segment, in the left side we use the first variable of the system Ox ei i 
and in the right side — the first variable of the system Ox ei ? . Consequently, for x = Oi t , we have 

^— e— (0) = — ^— e— (|ei|) . 
ox 1 ox 1 
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4.6 Justification: draft. 



Substituting the asymptotic solution of a form (|21[) we satisfy the convection-diffusion equation with a 
discrepancy O (e fc ) in L? norm. The boundary conditions are satisfied with the same accuracy. Then 
applying the Poincare-Friedrichs inequality for rod structures (see [I]) and the a priori estimate derived 
from the variational formulation, we get that 



where C does not depend on e. 

In particular the leading term cq satisfies the following estimate 



\\ c e - co|| L 2 (Be) < Ce. 



5 The partial asymptotic domain decomposition 

The constructed above asymptotic expansions of the solutions of the Stokes problem (fT]), (fj|, © and the 
diffusion-convection-sorption problem ([2]), ([5]), ([7]) allow us to apply the idea of the partial asymptotic 
domain decomposition [JJ UOj . We will cut off the two-dimensional subdomains of B £ containing the 
" stenosis areas" , bifurcations and eventually (for the Stokes problem) the entrance and exit elements 
by the lines orthogonal to the rectangles B e j at the distance S = Ke |lne| from the nodes (bifurcation 
points) and from the nodal points of the stenosis areas. We will call these subdomains the 2D zoom 
zones. Here K is independent of e and will be defined later. 

Then we pass to the ID description out of these subdomains. It means that we pass to the projection 
of the variational formulation of the Stokes problem on the Sobolev subspace of vector- valued functions 
having the Poiseuille "parabolic" shape out of these subdomains (see [HUH])- 

For the diffusion-convection-sorption problem we apply the projection on the Sobolev space of func- 
tions having vanishing derivatives of order greater than 2k in the direction orthogonal to B S j (also out 
of these 2D zoom zones). This choice of the projection space is motivated by Lemma [2] and formula 
(1141) . This gives us a variational formulation of the partially decomposed diffusion-convection problem 
and according to the general theory of the error estimate for the method of partial asymptotic domain 
decomposition ([HIE]]) we get the estimates: 

- for the difference of u £ and u de g (solution of the partially decomposed problem for Stokes equation) 
we get as in [8] that for any N there exists K independent of e such that, if 6 = Ke |lne| then 



<0(e N ) V~s; 



- for the difference of c e and c^f (solution of the partially decomposed problem for the diffusion- 
convection equation) we get that there exists K independent of e such that, if 8 = Ke |lne| then 

\\c e -c d j\\ H1{Be) <0{e k )V~e. (25) 
These estimates justify the application of the MAPDD. 

Remark 4 The interface conditions between the 2D parts of the domain and the ID parts follows from 
the variational formulation for the partially decomposed problem by integrating by parts. 

Remark 5 Although for k = estimate i25\) is not too precise, we will hold below a numerical experiment 
comparing the difference between the exact solution and the solution of partially decomposed problem 
in this simpliest case, when the projection space consists of functions with vanishing first transversal 
derivative out of the 2D zoom zones. 
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6 Numerical experiments 



Here we will compare the solution of the leading term ID equation (f22"|) - (j2"4|) to the numerical FEM 
solution of the coupled 2D flow-diffusion problem (fT|), ©, ©, ©, in a thin rectangle (0, 1) x (0, e). 
We will trace the boundary layer zones. 



O 1 

Fig. 4: Thin rectangle ("straight 
channel geometry") 



In the second part of this section we will discuss the numerical solution attained by the MAPDD for a 
one boundle tube structure corresponding to three segments e±, B2, e.%. 




Fig. 5: One bundle tube structun 



The finite clement discretization of the Stokes equations (Q} is based on the classical P2/P1 lagrangian 
finite element test functions in combination with the triangular finite element mesh. More precisely, the 
velocity field is approximated by quadratic lagrangian test functions while the pressure field is approx- 
imated with linear lagrangian test functions. As it is well known, this finite element flow formulation 
satisfies the classical Babuska-Brczzi compatibility condition and consequently produces numerically 
stable and adequate solution strategy for Stokes and Navier-Stokes problems. The concentration field c 
is approximated by quadratic lagrangian test functions. Finally, we would like to mention that for the 
studied flow conditions, no specific divergence problems were encountered. 



6.1 Straight channel geometry 

Before analyzing the case of a 2D bifurcation problem, we will present some results concerning the simple 
straight channel rectangle geometry B £ = (0, e) x (0, 1) . The main goal is to compare the predictions 
of our two methods (complete 2D and asymptotic ID model) in this simple flow conditions: e is taken 
equal to 0.05, viscosity /i = 1. 
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In this case, the velocity distribution is described by a planar Poiscuillc flow given by formula ([8|). 
The entry concentration is fixed as qo = 1, while the exit concentration is maintained to be q\ = 0.5. 
We calculate the solutions for (3 = 0.4 and three different values of the diffusion coefficient x. The 
comparison of 2D solution and the ID asymptotic solution is presented at Fig. 6. 

These variations of the diffusion coefficient are taken in order to find the limits of the asymptotic 
approximation. This asymptotic analysis was applied under the hypothesis of absence of other small 
parameters in the model. Indeed, when the diffusion coefficient becomes a second small parameter then 
the asymptotic analysis taking into account only one small parameter may be not too precise. 

D - QUI ; Beta ^ 11.4 : = 0,05 Mdx: 1,00 




-0.05 OJ05 0.1 0.15 02 025 03 035 0.4 0.45 0.5 055 0.6 0.65 0.7 0.75 OS 035 OS 035 1 1.05 



D * O.l ; Betd - 0.4 : eps — DH5 
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Mm:034» (l>)- 
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-0.05 0iJ5 0.1 0.15 02 025 03 035 0.4 0.45 0.5 0.55 0J5 065 0.7 0.75 OS 0J85 09 0S5 1 1J 

Mill: 0500 

Fig. 6: Distribution of the velocity field and distribution of the concentration for values 

(a) x = 0.01, (b) k = 0.1, (c) x = 1. 
Steady state concentration distribution: color map shows the concentration distribution 
the arrows show that the total flux density for the concentration distribution 
(i.e., diffusive plus convective flux). 



At low diffusion (Fig. 6(a) and 6(b)), the concentration distribution seems to be essentially 2D and 
therefore one can expect some differences between the corresponding ID and 2D approaches. On the 
other hand, when the diffusion increases, the concentration iso-levels become more and more planar 
and, consequently, we can expect better performance from the ID approach. This analysis is further 
confirmed by the direct comparison between the predictions of the ID and the 2D approaches presented 
in Fig. 7. 



17 



Concentration, c 






0=0.01 




0=0.1 




D=l 




D=10 




— -D=0X)1 




— 0=0.1 




* 0=1 




d=id 



Fig. 7: Direct comparison between the 2D (full line) and ID (dashed line) numerical 
predictions for 8 different values of the diffusion x. 
Concentration distribution: full line stands for the average over a cross-section of 
concentration distribution for the 2D geometry; 
dashed line and dots stand for the ID concentration distribution. 



As we can see from Fig. 7, some differences between ID and 2D predictions exist only for the value 
of x = 0.01 when the diffusion is 5 times smaller than e. So the numerical experiment confirms the great 
precision of the asymptotic solution (even in the case of small diffusion coefficient!). 



6.2 2D bifurcation geometry 

In this part, we extend our study to a more complex 2D bifurcation geometries. The 2D flow geometry 
is presented in Fig. 8(a). As it is seen from this figure, each channel could have different thickness and 
could be expressed in terms of so called streamline function ip according to the following definition: 



dip 
dx 



-u, 



ir 



dip 
dy 



U x 



which is calculated as a solution of the following differential equation 



Aip 



dU x dUy 



dy 



dx 



(26) 



(27) 



The flow kinematics around the bifurcation point is illustrated in Fig. 8(a). The corresponding pressure 
distribution is given in Fig. 8(b). As it is predicted by the asymptotic analysis, the pressure gradient in 
each arm is constant and naturally depends on the channel thickness and flow rate distribution. 
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Fig. 8: Model of the Stokes flow and the convection/diffusion process in a one boundle tube struc- 
ture: (a) flow kinematics, (b), (c) pressure distribution in the bifurcation zone, (d) and (e) iso- lines of 
concentration for two values of xr : x = 0.25 (d) and x = 0.10 (e), (f) comparison between the 2D and 
MAPDD solutions for the diffusion x = 0.1, 1 and 10. 
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It is important to emphasize that each geometry corner represents a singular point for the pressure, 
which can be seen at Fig. 8(c). These peaks are predicted by the corner singularities analysis f [41 [T3l 
[14j[T5]). Two typical solutions representing the cases of lower (x = 0.25) and higher (x = 10) diffusivity 
are given in Fig. 8(d) and 8(e). Like in the previously analyzed geometries, there is a critical diffusivity 
value, which ensures the validity of the simplified ID approach. In Fig. 8(f), we have given the direct 
comparison between ID and 2D predictions. As it is seen from these figures, the critical diffusivity value 
is around 1. 
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